/*
Copyright (c) 2013 Elizabeth Aston

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/


/* Diploid version of GeneticAlgorithm.cc, using initialisation type 3 (half on each peak with random pairs)
   Input arguments numGenotypes, mutationRate (X in 1000000), dominanceNum


- 	Two copies of the genome, one inherited from the mother, the other inherited from the father. Stored in a maternal array and a paternal array (where the sequence stored in position 1 of the maternal array is homologous to, 		and so pairs with, the sequence stored in position 1 of the paternal array). 
- 	Start with identical homologous sequences in each.

- 	Select 3 sequences at random from the population *.
- 	Mutate each individual (both constituent sequences) according to the per-base mutation rate.
- 	Work out the fitness of each of the resulting gametes **, and select which should be parents accordingly. There should be an equally small chance of either of the two fittest of the three being replaced (25%), and a larger 		chance of replacing the least fit (50%). 
- 	Crossover occurs between the the two constituent sequences of each individual. Randomly choose a cut point, and copy the maternal copy of the sequence up to the cut point, and the paternal copy from the cut point to the end. 		Store the resulting individual in a gamete array. The gamete array will now contain three strings representing three gametes produced by each of the three individuals.
-	At this stage there are now three sequences, each of which has undergone mutation and crossover, and two of which have been designated parents to replace the third.
-	Find the position corresponding to the third individual in the maternal and paternal arrays. Overwrite the sequence stored in the maternal array with the gamete of one of the parent individuals, and overwrite the sequence 		stored in the paternal array with the gamete of the other parent individual.

Go back to * and repeat until every individual has been selected once (and only once). The population has now gone through one generation.

DOMINANCE
1:	When calculating fitness, calculate the fitness of both the maternal and paternal sequences and use the HIGHER of the two as the fitness of the individual. This applies during reproduction (**) and when calculating the 		position of each individual in the fitness landscape.
2:	When calculating fitness, calculate the fitness of both the maternal and paternal sequences and use the LOWER of the two as the fitness of the individual. This applies during reproduction (**) and when calculating the 		position of each individual in the fitness landscape. 
3:	When calculating fitness, calculate the fitness of both the maternal and paternal sequences and use the (z*fM) + ((1-z)*fP) as the fitness (with z being a value >0 and <1) 
	This applies during reproduction (**) and when calculating the position of each individual in the fitness landscape.
*/


#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int numGenotypes;
int genotypeLength;
const int numGenerations = 10000;
const int numMutationsPerReproduction = 1;
int mutationRate; // probability (X in 1000000) of mutation
const int alphabetSize = 4; // Numerical representation of GATC
int **genotypesMaternal; /*[numGenotypes][genotypeLength]*/ // Stores the maternally inherited genes that make up the population
int **genotypesPaternal; /*[numGenotypes][genotypeLength]*/ // Stores the paternally inherited genes that make up the population
double *fitness;/*[numGenotypes]*/ // Stores the fitness of each genotype in the population
double fitness0 = 15;
double fitness1 = 10;
int radius0 = 2; // radius of X means up to X mutations away
int radius1 = 5;
int *target0; //[genotypeLength];
int *target1; //[genotypeLength];
int distanceBetweenTargets; //radius0+radius1 + 3; // keep in range [radius0+radius1,genotypeLength]
bool peak0;
bool peak1;
int **gameteArray;
int dominanceNum;
double z = 0.999999999999999; // Used with dominance type 3 to give a fitness of (z*fM) + ((1-z)*fP), where z >0 and <1

const int probabilityOfReplacing1stFittestOf3 = 250;
const int probabilityOfReplacing2ndFittestOf3 = 250;

int random(int n) 
{
  return (int) (((double) n)*rand()/(RAND_MAX+1.0));
}






void calculateIndividualFitness(int individual, bool checkOnPeak)
{
  /* Goes along the length of the current genotype
     Compares each bit with the corresponding position in both the target genotypes
     For every bit that matches the corresponding bit in the target, distanceToTarget is incremented
     Each distanceToTarget will therefore equal the distance away the current genotype is from each target genotype, in terms of bit differences */
  int maternalDistanceToTarget0 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypesMaternal[individual][i] != target0[i]) maternalDistanceToTarget0++;
  int maternalDistanceToTarget1 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypesMaternal[individual][i] != target1[i]) maternalDistanceToTarget1++;
  int paternalDistanceToTarget0 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypesPaternal[individual][i] != target0[i]) paternalDistanceToTarget0++;
  int paternalDistanceToTarget1 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypesPaternal[individual][i] != target1[i]) paternalDistanceToTarget1++;

  /* 
     DOMINANCE 1: Highest fitness.
     DOMINANCE 2: Lowest fitness.
     DOMINANCE 3: Average.
  */
  double maternalF0 = fitness0 * (1 - (maternalDistanceToTarget0/(double)radius0));
  double maternalF1 = fitness1 * (1 - (maternalDistanceToTarget1/(double)radius1));
  double paternalF0 = fitness0 * (1 - (paternalDistanceToTarget0/(double)radius0));
  double paternalF1 = fitness1 * (1 - (paternalDistanceToTarget1/(double)radius1));


  /* 
    maternalF0 is the fitness of the maternal sequence based on peak 0, and maternalF1 is the fitness of the maternal sequence based on peak 1.
    The higher the fitness value, the closer the individual is to the top of a peak. If maternalF0 is higher than maternalF1, the maternal sequence of the 
    individual is closer to peak 0 (and vice versa). Take the highest of the two values and store it in mF. This is the fitness of the maternal sequence.
    Do the same for the paternal fitnesses based on peak 0 and 1, to give a final paternal fitness pF.
    N.B. The minimum fitness is 0, therefore if either mF or pF have a negative value, set them to 0
  */

  double mF = maternalF0;
  
  if(maternalF1 > maternalF0)
  {
    mF = maternalF1;
  } 
  
  if(mF < 0)
  {
    mF = 0;
  }  


  double pF = paternalF0;
  
  if(paternalF1 > paternalF0)
  {
    pF = paternalF1;
  } 

  if(pF < 0)
  {
    pF = 0;
  }

  /*
    There is now a fitness value for the maternal sequence (mF), and a fitness value for the paternal sequence (pF).
    Use one of dominance mechanisms 1, 2 or 3 to select a final overall fitness value for this individual.
  */

  // Choose the sequence with the highest fitness
  if(dominanceNum == 1)
  {
    if(mF > pF)
    {
      fitness[individual] = mF;
    }
    else
    {
      fitness[individual] = pF;
    }
  }

  // Choose the sequence with the lowest fitness
  if(dominanceNum == 2)
  {
    if(mF < pF)
    {
      fitness[individual] = mF;
    }
    else
    {
      fitness[individual] = pF;
    }
  }

  // Choose a fitness somewhere in between (determined by the value of z (where z is >0 and <1)
  // N.B. If z = 0.5, it will be the average fitness of mF and pF
  if(dominanceNum == 3)
  {
    // Low values of z should lean towards the minimum fitness, while high values of z should be closer to the maximum fitness
    // First work out whether mF or pF have the highest fitness

    double max = pF;
    double min = mF;

    if(mF > pF)
    {
      max = mF;
      min = pF;
    }

    fitness[individual] = (z*max) + ((1-z)*min);
  }





  if(checkOnPeak == true)
  {
    if((maternalDistanceToTarget0 < radius0) || (paternalDistanceToTarget0 < radius0)) 
    {
      peak0 = true;
    }
    if((maternalDistanceToTarget1 < radius1) || (paternalDistanceToTarget1 < radius1)) 
    {
      peak1 = true;
    }
  }
}





void initialise()
{
  unsigned int s; ifstream u("/dev/urandom",ios::in|ios::binary);
u.read((char*)&s,sizeof(s)); if (!u) {cerr << "Failed to get seed from kernel random number source device" << endl; exit(1);}
  srand(s);

  /* Initialise target0 and target1 (which represent the two peaks)
     target0 stores the genotype at the top of peak 0
     target1 stores the genotype at the top of peak 1 */
  for(int i = 0; i < genotypeLength; i++)
  {
    target0[i] = 0; // target0 is a sequence of 0s
    target1[i] = target0[i]; // target1 is initialised to all 0s but will be mutated according to the value of distanceBetweenTargets in the next step
  }


  /* The following code mutates target1 a number of times equal to distanceBetweenTargets
     This ensures target1 is the correct number of bits different in fitness to target0 */

  /* change stores true or false for each bit, to ensure that once a bit gets mutated it cannot be mutated back
     i.e., preventing back mutation */
  bool change[genotypeLength]; 

  // Initialises each position in the change array to false
  for(int i = 0; i < genotypeLength; i++) 
  {
    change[i] = false;
  }

  /* Changes a number of bits in target1 to make it the required number of bits
     different to target0 */
  for(int d = 0; d < distanceBetweenTargets; d++) 
  {
    int i;
    
    do 
    {
      i = random(genotypeLength); // Randomly chooses a position in the genotype that hasn't already been changed
    } 
    while(change[i] != false);

    change[i] = true; // Sets change to true to ensure this position cannot be selected again
    
    int base;
    
    do
    {
      base = random(alphabetSize); // Randomly generates a new value to go in the selected position of the genotype
    } 
    while(base == target1[i]); // Ensures the new value isn't the same as the one already stored in that position of the target1 genotype

    target1[i] = base; // Replaces the existing value with the mutant
  }



  // Generate the initial population
  for(int i = 0; i < numGenotypes; i++)
  {
    for(int j = 0; j < genotypeLength; j++)
    {
      genotypesMaternal[i][j] = target0[j]; // Initialises the whole population to target0
      genotypesPaternal[i][j] = target0[j]; // Initialises the whole population to target0
    }
  }

  bool *chosen;
  chosen = new bool[numGenotypes];

  // Initialise chosen to false
  for(int i = 0; i < numGenotypes; i++)
  {
    chosen[i] = false;
  }

  // Puts half of the maternal sequences on to peak 1
  for(int i = 0; i < numGenotypes/2; i++)
  {
    int randomIndividual;

    // Picks an individual at random that hasn't been previously chosen
    do
    {
      randomIndividual = random(numGenotypes);
    }
    while(chosen[randomIndividual]);

    chosen[randomIndividual] = true;

    // Places this individual on peak 1
    for(int j = 0; j < genotypeLength; j++)
    {
      genotypesMaternal[randomIndividual][j] = target1[j]; 
    }
  }

  // Reinitialises chosen array
  for(int i = 0; i < numGenotypes; i++)
  {
    chosen[i] = false;
  }

  // Puts half of the paternal sequences on to peak 1
  for(int i = 0; i < numGenotypes/2; i++)
  {
    int randomIndividual;

    // Picks an individual at random that hasn't been previously chosen
    do
    {
      randomIndividual = random(numGenotypes);
    }
    while(chosen[randomIndividual]);

    chosen[randomIndividual] = true;

    // Places this individual on peak 1
    for(int j = 0; j < genotypeLength; j++)
    {
      genotypesPaternal[randomIndividual][j] = target1[j]; 
    }
  }

  // Calculates the updated fitness of each individual
  for(int j = 0; j < numGenotypes; j++) 
  {
    calculateIndividualFitness(j, false);
  }

  for(int i = 0; i < 2; i++)
  {
    for(int j = 0; j < genotypeLength; j++)
    { 
      gameteArray[i][j] = 0;
    }
  }

  delete[] chosen;
}




/* Crossover occurs between the maternal and paternal sequences of an individual
   The resulting sequence is stored in the gameteArray */
void crossover(int individual, int gamete)
{
  // Randomly generates a position in the genotype at which to cut (0 < cutPoint > genotypeLength to ensure crossover occurs)
  int cutPoint = 1 + random(genotypeLength - 1);

  for(int g = 0; g < cutPoint; g++)
  {
    gameteArray[gamete][g] = genotypesMaternal[individual][g]; // Copies the maternally inherited sequence of the current individual up to the cut point
  }  
  for(int g = cutPoint; g < genotypeLength; g++)
  {
    gameteArray[gamete][g] = genotypesPaternal[individual][g]; // Copies the maternally inherited sequence of the current individual from the cut point to the end of the sequence
  }
}




/* Mutate each position in the individual sequence with a given probability (mutationRate)
   There is an equal chance of mutation to any of the other possible values (GATC) */
void mutate(int gamete)
{
  int r = 0; 
  int newValue;

  /* Goes along each position of the sequence in the gamete of the current individual */
  for(int j = 0; j < genotypeLength; j++)
  {
    r = random(1000000);

    // If the mutation rate is greater than or equal to random value r, mutate the current position
    if(mutationRate > r)
    {
      /* Do the point mutation by randomly generating a new value to replace the original value stored in the gametes array
         N.B. The new value cannot be the same as the original value (otherwise there has been no mutation) */

      newValue = gameteArray[gamete][j]; // Copies the original value into newValue, which is then mutated

      while(newValue == gameteArray[gamete][j]) // Makes sure the mutation is done again if necessary
      {
        newValue = random(alphabetSize);
      }

      gameteArray[gamete][j] = newValue; // Replaces the original value in the gametes array with the new mutant value
    }
  }
}




void nextStep()
{
  int *popArray;
  popArray = new int[numGenotypes];

  // Initialises popArray as an ordered list of 0 up to numGenotypes-1
  for(int i = 0; i < numGenotypes; i++)
  {
    popArray[i] = i;
  }

  // Sets the initial position to be at the top of the list
  int position = numGenotypes - 1;

  // Goes along the length of popArray and changes the order to random
  for(int i = 0; i < numGenotypes; i++)
  {
    int replace = random(position+1); // Generates a random number between 0 and the current position

    // Switches the value stored in position with the value at the randomly selected position
    int temp = popArray[position];
    popArray[position] = popArray[replace];
    popArray[replace] = temp;

    // Moves position down one in the list
    position--;
  }


  int a,b,c;

  for(int i = 0; i < numGenotypes-2; i+=3)
  {
    a = popArray[i];
    b = popArray[i+1];
    c = popArray[i+2]; 
    

    // Reorder such that c is the least fit
    if(fitness[a] < fitness[c]) 
    {
      int temp = a;
      a = c;
      c = temp; 
    }
    if(fitness[b] < fitness[c]) 
    { 
      int temp = b;
      b = c;
      c = temp;
    }

    // Make sure a is the fittest
    if(fitness[a] < fitness[b])
    {
      int temp = a;
      a = b;
      b = temp;
    }

    int child;
#define  nonStandardSteadyStageGA
#ifndef  nonStandardSteadyStageGA
    child = c;
#else  //nonStandardSteadyStageGA
  #warning "Non-standard Steady Stage GA: individual to be replaced chosen by rank based selection where probabilities of selecting fittest and second fittest of 3 may be non-zero"
    int r = random(1000);
    if (r<probabilityOfReplacing1stFittestOf3) {child = a; a = c;} // Copy c into a?
    else if
    (r<(probabilityOfReplacing1stFittestOf3+probabilityOfReplacing2ndFittestOf3)) {child = b; b = c;} // Copy c into b?
    else child = c;
    // nb may well be that child==a or child==b but standard crossover code (below and in crossover()) can cope with that
#endif

    /* Crossover the constituent maternal and paternal sequences of individual a and b
       The parameters are individual and gamete, i.e., crossover of individual a produces the gamete stored in position 0 of the gamete array */
    int gameteA = 0;
    int gameteB = 1;

    crossover(a, gameteA);
    crossover(b, gameteB);

   // Mutate the resulting gametes. This simulates the errors due to DNA replication.
    for(int m = 0; m < numMutationsPerReproduction; m++) 
    {
      mutate(gameteA);
      mutate(gameteB);
    }     


    /* At this stage there are now two gametes, each of which has undergone mutation and crossover, and an individual which has been designated the child. 
       Find the position corresponding to the child individual in the maternal and paternal arrays. Overwrite the sequence stored in the maternal array with the gamete 
       of one of the parent individuals, and overwrite the sequence stored in the paternal array with the gamete of 
       the other parent individual.
    */
    if(random(2) == 1) // Randomly generate either a 2 or a 1. If it's 1, switch a and b. This ensures it is random which of the parent gametes is designated maternal and which is paternal             
    {
      int tempG = gameteA;
      gameteA = gameteB;
      gameteB = tempG;
    }

    /* Overwrite the child individual with the gamete produced by each parent
       GameteA is designated the maternal sequence, and gameteB is designated the paternal sequence */
    for(int g = 0; g < genotypeLength; g++)
    {
      genotypesMaternal[child][g] = gameteArray[gameteA][g];
      genotypesPaternal[child][g] = gameteArray[gameteB][g];
    }  

    // Recalculate the fitness of the child individual (false stops onPeak being called unnecessarily)
    calculateIndividualFitness(child, false);  
  }

  delete[] popArray;
}





int main(int argc, char* argv[])
{ 
  #warning "Input arguments are numGenotypes, genotypeLength, distanceBetweenTargets, mutationRate (X in 1000000), dominanceNum (1 highest fitness, 2 lowest fitness, 3 average fitness)"

  if(argv[1]==0) {cerr << "Missing argument" << endl; exit(1);}
  numGenotypes = atoi(argv[1]);
  
  fitness = new double[numGenotypes];

  if(argv[2]==0) {cerr << "Missing argument" << endl; exit(1);}
  genotypeLength = atoi(argv[2]);

  target0 = new int[genotypeLength];
  target1 = new int[genotypeLength];

 genotypesMaternal = new int*[numGenotypes];
  for (int g=0;g<numGenotypes;g++)
     genotypesMaternal[g] = new int[genotypeLength];
  genotypesPaternal = new int*[numGenotypes];
  for (int g=0;g<numGenotypes;g++)
     genotypesPaternal[g] = new int[genotypeLength];

  // Creates an array to store the gametes produced by crossover of each individual
     gameteArray = new int*[2];
     for(int g = 0; g < 2; g++){gameteArray[g] = new int[genotypeLength];}

  if(argv[3]==0) {cerr << "Missing argument" << endl; exit(1);}
  distanceBetweenTargets = atoi(argv[3]);

  if(argv[4]==0) {cerr << "Missing argument" << endl; exit(1);}
  mutationRate = atoi(argv[4]);

  if(argv[5]==0) {cerr << "Missing argument" << endl; exit(1);}
  dominanceNum = atoi(argv[5]);
  if((dominanceNum > 3) || (dominanceNum < 1)) {cerr << "Dominance options 1, 2 or 3 only" << endl; exit(1);}


  initialise();

  int peak0LostAt = -1;
  int peak1LostAt = -1;

  //double start_time = clock();

  for(int generation = 1; generation <= numGenerations; generation++)
  {
    peak0 = peak1 = false;
    nextStep();

    // Calculates the fitness of each individual in the population and works out whether there are any individuals on either peak
    for(int i = 0; i < numGenotypes; i++)
    {
      calculateIndividualFitness(i, true); // True tells calculateIndividualFitness to call the onPeak method
    }

   
    // If it's the first time the peak has been lost, record the current generation
    if ((peak0==false) && (peak0LostAt==-1)) peak0LostAt = generation;
    if ((peak1==false) && (peak1LostAt==-1)) peak1LostAt = generation; 
  
    // If both peaks have been lost, output when they were lost and end the program
    if((peak0LostAt > -1) && (peak1LostAt > -1)) break; 
  }

  cout << peak0LostAt << "\t" << peak1LostAt << endl;


  delete[] fitness;
  delete[] target0;
  delete[] target1;
  for (int g=0;g<numGenotypes;g++)
     delete[] genotypesMaternal[g];
  delete[] genotypesMaternal;
  for (int g=0;g<numGenotypes;g++)
     delete[] genotypesPaternal[g];
  delete[] genotypesPaternal;
  for(int g = 0; g < 3;g++)
    delete[] gameteArray[g];
  delete[] gameteArray;
  
  return 0;
}

